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ABSTRACT 

We simulate the growth of galaxies and their central supermassive black holes by 
implementing a suite of semi-analytic models on the output of the Millennium Run, 
a very large simulation of the concordance ACDM cosmogony. Our procedures fol- 
low the detailed assembly history of each object and are able to track the evolution 
of all galaxies more massive than the Small Magellanic Cloud throughout a volume 
comparable to that of large modern redshift surveys. In this first paper we supple- 
ment previous treatments of the growth and activity of central black holes with a new 
model for ‘radio’ feedback from those AGN that lie at the centre of a quasistatic X-ray 
emitting atmosphere in a galaxy group or cluster. We show that for energetically and 
observationally plausible parameters such a model can simultaneously explain: (1) the 
low observed mass drop-out rate in cooling flows; (ii) the exponential cut-off at the 
bright end of the galaxy luminosity function; and (iii) the fact that the most massive 
galaxies tend to be bulge-dominated systems in clusters and to contain systematically 
older stars than lower mass galaxies. This success occurs because static hot atmo- 
spheres form only in the most massive structures, and radio feedback (in contrast, for 
example, to supernova or starburst feedback) can suppress further cooling and star 
formation without itself requiring star formation. We discuss possible physical mod- 
els which might explain the accretion rate scalings required for our phenomenological 
‘radio mode’ model to be successful. 
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scale structure (Peacock et al/|2001), the baryon fraction in 
rich clusters (White et al|[1993) and the theory of Big Bang 
nucleosynthesis (Olive et a1/[2000). A working model for the 
growth of all structure thus appears well established. 


1 INTRODUCTION 


The remarkable agreement between recent measurements of 
cosmic structure over a wide range of length- and time- 


scales has established a standard paradigm for structure 
formation, the ACDM cosmogony. This model can simulta- 
neously match the microwave background fluctuations seen 
at z ~ 1000 (e.g. Breed albo the power spectrum 
of the low redshift galaxy distribution (e.g. [Percival et ail 
[2003 2004), the nonlinear mass distribu- 
tion at low redshift as characterised by cosmic shear (e.g. 
2002) and the structure seen in the 
z —3 Ly a forest (e.g. MEL the ste 2003). It also re- 
produces the present acceleration of the cosmic expansion in- 
ferred from supernova observations (Perlmutter et al 
[Riess ot arTrosgh, and it is consistent with the mass budget 


inferred for the present universe from the dynamics of large- 


In this cosmogony, galaxies form when gas condenses at 
the centres of a hierarchically merging population of dark 
haloes, as originally proposed bylWhin E Reed (1978). At- 
tempts to understand this process in detail have consistently 
run into problems stemming from a mismatch in shape be- 
tween the predicted distribution of dark halo masses and 
the observed distribution of galaxy luminosities. Most stars 
are in galaxies of Milky Way brightness; the galaxy abun- 
dance declines exponentially at brighter luminosities and in- 
creases sufficiently slowly at fainter luminosities that rela- 
tively few stars are in dwarfs. In contrast, the theory pre- 
dicts a much broader halo mass distribution — a constant 
mass-to-light ratio would produce more high and low lu- 
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minosity galaxies than are observed while underpredicting 
the number of galaxies like the Milky Way. Attempts to 
solve these problems initially invoked cooling inefficiencies 
to reduce gas condensation in massive systems, and super- 
nova feedback to reduce star formation efficiency in low mass 
systems (White & Reed1978 1991). Forma- 
tion of dwarfs may also be suppressed by photoionisation 
heating (Efstathioul1992). As|Thoul & Weinberg em- 
phasised, cooling effects alone are too weak to produce the 
bright end cut-off of the luminosity function, and recent at- 
tempts to fit observed luminosity functions have been forced 
to include additional feedback processes in massive systems 
(e.g. [Benson et al][2003). In this paper we argue that ra- 
dio sources may provide the required feedback while at the 
same time providing a solution to two other long-standing 
puzzles. 

An important unanswered question is why the gas at 
the centre of most galaxy clusters is apparently not con- 
densing and turning into stars when the observed X-ray 
emission implies a cooling time which is much shorter than 
the age of the system. This cooling flow puzzle was noted 
as soon as the first X-ray maps of clusters became avail- 
able (Cowie & Binacd[i972 a 
was made more acute when X-ray spectroscopy demon- 
strated that very little gas is cooling through temperatures 
even a factor of three below that of the bulk of the gas 
(Peterson et all[2001} [Tamura et al][2001].. A clue to the so- 
lution may come from the observation (Burns et al] 
that every cluster with a strong cooling flow also contains 
a massive and active central radio galaxy. 
suggested that radio galaxies might regulate cool- 
ing flows, and this idea has gained considerable recent 
support from X-ray maps which show direct evidence for 
an interaction between radio lobes and the intracluster 


gas (Fabian et al] [2003: [McNamara et al] [200d, [2003). A 


number of authors have suggested ways in which the ra- 


dio source might replace the thermal energy lost to X- 
ray emission (Binney & Tabo ; |Churazov et al. 


Brüggen & Kaiser |2002; [Ruszkowski & Begelman 
; . We do not focus 
on this aspect of the problem here, but rather demonstrate 
that if the scaling properties of radio source feedback are set 
so they can plausibly solve the cooling flow problem they 
induce a cut-off at the bright end of the galaxy luminosity 
function which agrees well with observation. 


Another puzzling aspect of the galaxy population is 
the fact that the most massive galaxies, typically ellipti- 
cals in clusters, are made of the oldest stars and so fin- 
ished their star formation earlier than lower mass galax- 
ies (Bender & Saglid [1999). Confirming evidence for this 
comes from look-back studies which show that both star- 
formation and AGN activity take place more vigorously 
and in higher mass objects at redshifts of 1 to 2 than in 


the present Universe e.g. [Shaver et all[199&: [Madau et al] 
Lo [Cowie et all (1996) termed this phenomenon ‘down- 
sizing’, and prima facie it conflicts with hierarchical growth 
of structure in a ACDM cosmogony where massive dark 
haloes assemble at lower redshift than lower mass haloes 
(e.g. [Lacey & Cold[1993). 'This puzzle is related to the pre- 
vious two; the late-forming high mass haloes in ACDM corre- 
spond to groups and clusters of galaxies, and simple theories 
predict that at late times their central galaxies should grow 


to masses larger than those observed through accretion from 
cooling flows. In the model we present below, radio galaxies 
prevent significant accretion, thus limiting the mass of the 
central galaxies and preventing them from forming stars at 
late times when their mass and morphology can still change 
through mergers. The result is a galaxy luminosity function 
with a sharper high-mass cut-off in which the most massive 
systems are red, dead and elliptical. 

To make quantitative predictions for the galaxy popu- 
lation in a ACDM universe it is necessary to carry out sim- 
ulations. Present numerical capabilities allow reliable sim- 
ulation of the coupled nonlinear evolution of dark mat- 
ter and diffuse gas, at least on the scales which determine 
the global properties of galaxies. Once gas cools and con- 
denses into halo cores, however, both its structure and the 
rates at which it turns into stars and feeds black holes are 
determined by small-scale ‘interstellar medium’ processes 
which are not resolved. These are usually treated through 
semi-analytic recipes, parameterised formulae which encap- 
sulate ‘subgrid’ physics in terms of star formation thresh- 
olds, Schmidt ‘laws’ for star formation, Bondi models for 
black hole feeding, etc. The form and the parameters of 
these recipes are chosen to reproduce the observed sys- 
tematics of star formation and AGN activity in galaxies 
(e.g. 1998). With a well-constructed scheme it 
is possible to produce stable and numerically converged 
simulations which mimic real star-forming galaxies remark- 


ably well (Springel & Hernquist |2003a)). In strongly star- 


forming galaxies, massive stars and supernovae produce 
winds which redistribute energy, mass and heavy elements 
over large regions [1990; [1999). 
Even stronger feedback is possible, in principle, from AGN 
(Beschman er ai]frognp. In both cases the determining pro- 
cesses occur on very small scales and so have to be included 
in simulations through parametrised semi-analytic models. 
Unfortunately, the properties of simulated galaxies turn out 
to depend strongly on how these unresolved star-formation 
and feedback processes are treated. 

Since the diffuse gas distribution and its cooling onto 
galaxies are strongly affected by the description adopted for 
the subgrid physics, every modification of a semi-analytic 
model (or of its parameters) requires a simulation to be re- 
peated. This makes parameter studies or tests of, say, the ef- 
fects of different AGN feedback models into a very expensive 
computing exercise. A cost-effective alternative is to repre- 
sent the behaviour of the diffuse gas also by a semi-analytic 
recipe. Since the dark matter couples to the baryons only 
through gravity, its distribution on scales of galaxy haloes 
and above is only weakly affected by the details of galaxy for- 
mation. Its evolution can therefore be simulated once, and 
the evolution of the baryonic component can be included 
in post-processing by applying semi-analytic models to the 
stored histories of all dark matter objects. Since the sec- 
ond step is computationally cheap, available resources can 
be used to carry out the best possible dark matter simula- 
tion, and then many parameter studies or tests of alterna- 
tive models can be carried out in post-processing. This sim- 
ulation approach was first introduced by kaut aan ct all 
and it is the approach we apply here to the Mil- 
lennium Run, the largest calculation to date of the evo- 
lution of structure in the concordance ACDM cosmogony 


(Springel et alJ[2005b)). 
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This paper is organised as follows. In Section B] we de- 
scribe the Millennium Run and the post-processing we car- 
ried out to construct merging history trees for all the dark 
haloes within it. SectionB]presents the model for the forma- 
tion and evolution of galaxies and their central supermas- 
sive black holes that we implement on these merging trees. 
Section Ø] describes the main results of our modelling, con- 
centrating on the influence of ‘radio mode’ feedback on the 
properties of the massive galaxy population. In Section [5] 
we discuss physical models for black hole accretion which 
may explain the phenomenology required for our model to 
be successful. Finally, Section [6]summarises our conclusions 
and suggests some possible directions for future investiga- 
tion. 


2 THE DARK MATTER SKELETON: THE 
MILLENNIUM RUN 


Our model for the formation and evolution of galaxies and 
their central supermassive black holes is implemented on 
top of the Millennium Run, a very large dark matter simu- 
lation of the concordance ACDM cosmology with 2160? ~ 
1.0078 x 10° particles in a periodic box of 500 ^^! Mpc on 
a side. A full description is given in[Springel et al] (2005); 
here we summarise the main simulation characteristics and 
the definition and construction of the dark matter merging 
history trees we use in our galaxy formation modelling. The 
dark matter distribution is illustrated in the top panel of 
Fig. [I] for a 330 x 280 x 15 ^! Mpc slice cut from the full 
volume. The projection is colour coded by density and lo- 
cal velocity dispersion, and illustrates the richness of dark 
matter structure for comparison with structure in the light 
distribution to which we will come later. Dark matter plots 
on a wider range of scales may be found in Bornedal] 
(20050). 


2.1 Simulation characteristics 


We adopt cosmological parameter values consistent with 
a combined analysis of ty eas dlls et all 
and first-year WMAP data (Spergel _et_al./2003;|Seljak et al 
[2005). They are Qm = Qam + Q» = 0.25, Qp = 0.045, 
h = 0.73, Qa = 0.75, n = 1, and og = 0.9. Here Qm de- 
notes the total matter density in units of the critical density 
for closure, perit = 3H /(81G). Similarly, Q» and Q4 denote 
the densities of baryons and dark energy at the present day. 
The Hubble constant is given as Ho = 100 kms"! Mpc "^ !, 
while og is the rms linear mass fluctuation within a sphere 
of radius 8h~'Mpc extrapolated to z —0. 

The chosen simulation volume is a periodic box of 
size 500h~'Mpc, which implies a particle mass of 8.6 x 
105 A^! Ma. This volume is large enough to include inter- 
esting objects of low space density, such as quasars or rich 
galaxy clusters, the largest of which contain about 3 mil- 
lion simulation particles at z — 0. At the same time, the 
mass resolution is sufficient that haloes that host galaxies 
as faint as 0.1 L, are typically resolved with at least ~ 100 
particles. Note that although discreteness noise significantly 
affects the merger histories of such low mass objects, the 
galaxies that reside in halos with £100 particles are usu- 
ally sufficiently far down the luminosity function that any 


uncertainty in their properties has little impact on our re- 
sults or conclusions. 

The initial conditions at z — 127 were created by dis- 
placing particles from a homogeneous, ‘glass-like’ distribu- 
tion using a Gaussian random field with a 
ACDM linear power spectrum as given by the Boltzmann 
code CMBFAST Elis & zadar [1999). The simula- 
tion was then evolved to the present epoch using a leapfrog 
integration scheme with individual and adaptive time steps, 
with up to 11 000 time steps for individual particles. 

The simulation itself was carried out with a spe- 
cial version of the GADGET-2 code 
optimised for very low memory consump- 
tion so that it would fit into the nearly 1 TB of physi- 


cally distributed memory available on the parallel IBM p690 
computer! used for the calculation. The computational al- 
1995; 


gorithm used the ‘TreePM’ method 
2000 2002) to evaluate gravitational forces, combin- 


ng a hierarchical multipole expansion, or ‘tree’ algorithm 
(bano: Hud i, and a classical, Fourier transform 
particle-mesh method 1981). An ex- 
plicit force-split in Fourier-space produces a very nearly 
isotropic force law with negligible force errors at the force 
matching scale. The short-range gravitational force law is 
softened on comoving scale 5 ^! kpc (Plummer-equivalent) 
which may be taken as the spatial resolution limit of the 
calculation, thus achieving a dynamic range of 10? in 3D. 
'The calculation, performed in parallel on 512 processors, 
required slightly less than 350000 processor hours of CPU 
time, or 28 days of wall-clock time. 


2.2 Haloes, substructure, and merger tree 
construction 


Our primary application of the Millennium Run in this pa- 
per uses finely resolved hierarchical merging trees which en- 
code the full formation history of tens of millions of haloes 
and the subhaloes that survive within them. These merging 
history trees are the backbone of the model that we imple- 
ment in post-processing in order to simulate the wide range 
of baryonic processes that are important during the forma- 
tion and evolution of galaxies and their central supermassive 
black holes. 

We store the full particle data between z —20 and z=0 
at 60 output times spaced in expansion factor according to 
the formula 


log (1 + zn) = n (n + 35)/4200 . (1) 


This spacing is ‘locally’ logarithmic but becomes smoothly 
finer at lower redshift, with a temporal resolution by red- 
shift zero of approximately 300 Myr. Additional outputs are 
added at z — 30, 50, 80, 127 to produce a total of 64 snap- 
shots in all. We note that each snapshot has a total size in 
excess of 300 GB, giving a raw data volume of nearly 20 TB. 

'Together with each particle coordinate dump, the sim- 
ulation code directly produces a friends-of-friends (FOF) 
group catalogue on the fly and in parallel. FOF groups are 


1 This computer is operated by the Computing Centre of the 
Max-Planck Society in Garching, Germany. 
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125 Mpo/h 


Figure 1. The redshift zero distribution of dark matter (top) and of galaxy light (bottom) for a slice of thickness 15 ^ -  Mpc, cut from 
the Millennium Run. For the dark matter distribution, intensity encodes surface density and colour encodes local velocity dispersion. 
For the light distribution, intensity encodes surface brightness and colour encodes mean B—V colour. The linear scale is shown by the 
bar in the top panel. 
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defined as equivalence classes where any pair of two parti- 
cles is placed into the same group if their mutual separation 
is less than 0.2 of the mean particle separation 
[1983). 'This criterion combines particles into groups with 
a mean overdensity of about 200, corresponding approxi- 
mately to that expected for a virialised group. The group 
catalogue saved to disk for each output only kept groups 
with at least 20 particles. 

High-resolution simulations like the present one ex- 
hibit a rich substructure of gravitationally bound dark mat- 
ter subhaloes orbiting within larger virialised halos (e.g. 
[Ghigna et al] [1998). The FOF group-finder built into our 
simulation code identifies the haloes but not their subha- 
los. Since we wish to follow the fate of infalling galaxies and 
their halos, and these are typically identifiable for a sub- 
stantial time as a dark matter subhalo within a FOF halo, 
we apply in post-processing an improved and extended ver- 
sion of the SUBFIND algorithm of [Springel et al] (20014). 
This computes an adaptively smoothed dark matter density 
field within each halo using a kernel-interpolation technique, 
and then exploits the topological connectivity of excursion 
sets above a density threshold to identify substructure can- 
didates. Each substructure candidate is subjected to a grav- 
itational unbinding procedure. If the remaining bound part 
has more than 20 particles, the subhalo is kept for further 
analysis and some of its basic physical properties are deter- 
mined (angular momentum, maximum of its rotation curve, 
velocity dispersion, etc.). After all subhaloes are identified 
they are extracted from the FOF halo so that the remain- 
ing featureless ‘background’ halo can also be subjected to 
the unbinding procedure. This technique, however, neglects 
the fact that substructures embedded within a halo help to 
bind its material, and thus removing them can, in princi- 
pal, unbind some of the FOF halo that may otherwise still 
be loosely bound. We accept this small effect for techni- 
cal reasons related to the robustness of our halo definition 
procedures and the need for unambigous particle-subhalo 
assignments in our data structures. We note that the to- 
tal mass in substructures is typically below 1096 and often 
substantially smaller. Thus any bias in the bound mass of 
the parent halo due to additional unbinding is always very 
small. 

'To compute a virial mass estimate for each FOF halo we 
use the spherical-overdensity approach, where the centre is 
determined using the minimum of the gravitational potential 
within the group and we define the boundary at the radius 
which encloses a mean overdensity of 200 times the critical 
value. The virial mass, radius and circular velocity of a halo 
at redshift z are then related by 


3 
Wir 
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Mir = H?(z) BR = (2) 
G 
where H(z) is the Hubble constant at redshift z. 

At z=0 our procedures identify 17.7 x 10° FOF groups, 
down from a maximum of 19.8 x 10° at z — 1.4 when groups 
were more abundant but of lower mass on average. The z —0 
groups contain a total of 18.2 x 10° subhaloes, with the 
largest FOF group containing 2328 of them. (Note that with 
our definitions, all FOF groups contain at least one subhalo, 
the main subhalo which is left over after removal of any 
substructure and the unbound component. We require this 
main subhalo to contain at least 20 particles.) 


Having found all haloes and subhaloes at all output 
snapshots, we then characterise structural evolution by 
building merging trees that describe in detail how haloes 
grow as the universe evolves. Because structures merge hier- 
archically in CDM universes there can be several progenitors 
for any given halo, but in general there is only one descen- 
dant. (Typically the cores of virialised dark matter struc- 
tures do not split into two or more objects.) We therefore 
construct merger trees based on defining a unique descen- 
dant for each halo and subhalo. This is, in fact, sufficient to 
define the entire merger tree, since the progenitor informa- 
tion then follows implicitly. Further details can be found in 
Springel et al] (2005b). 

We store the resulting merging histories tree by tree. 
Since each tree contains the full formation history of some 
z=0 halo, the physical model for galaxy formation can be 
computed sequentially tree by tree. It is thus unnecessary 
to load all the history information on the simulation into 
computer memory at the same time. Actually, this would 
be currently impossible, since the combined trees contain a 
total of around 800 million subhaloes for each of which a 
number of attributes are stored. Note that, although evolv- 
ing the galaxy population sequentially in this way allows us 
to consider, in principle, all interactions between galaxies 
that end up in the same present-day FOF halo, it does not 
allow us to model longer range interactions which might take 
place between galaxies that end up in different FOF halos 
at z=0. 


3 BUILDING GALAXIES: THE 
SEMI-ANALY TIC MODEL 


3.1 Overview 


In the following sub-sections we describe the baryonic 
physics of one particular model for the formation and evo- 
lution of galaxies and of their central supermassive black 
holes. A major advantage of our simulation strategy is that 
the effects of parameter variations within this model (or 
indeed alternative assumptions for some of the processes) 
can be explored at relatively little computational expense 
since the model operates on the stored database of merger 
trees; the simulation itself and the earlier stages of post- 
processing do not need to be repeated. We have, in fact, 
explored a wide model and parameter space to identify our 
current best model. We summarise the main parameters of 
this model, their ‘best’ values, and their plausible ranges in 
Table [I] These choices produce a galaxy population which 
matches quite closely many observed quantities. In this pa- 
per we discuss the field galaxy luminosity-colour distribu- 
tion; the mean stellar mass — stellar age relation; the Tully- 
Fisher relation, cold gas fractions and gas-phase metallicities 
of Sb/c spirals; the colour — magnitude relation of ellipticals; 
the bulge mass — black hole mass relation; and the volume- 
averaged cosmic star-formation and black hole accretion his- 
tories. In[Springel et al] we also presented results for 
galaxy correlations as a function of absolute magnitude and 
colour, for the baryonic ‘wiggles’ in the large-scale power 
spectrum of galaxies, and for the abundance, origin and 
fate of high redshift supermassive black holes which might 


correspond to the z ~ 6 quasars discovered by the SDSS 
(Enn et allogi) 
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Table 1. A summary of our main model parameters and their best values and plausible ranges, as described in the text. Once set, these 
values are kept fixed for all results presented in this paper, in particular for models in which AGN feedback is switched off. 


parameter description best value plausible range 
fb cosmic baryon fraction (3:3) 0.17 fixed 
20, Zr redshift of reionization ({3.3) 8,7 fixed 
fBH merger cold gas BH accretion fraction (B41) 0.03 0.02 — 0.04 
KAGN quiescent hot gas BH accretion rate (Mo yr-!) (B43) 6x10-6  (4—8) x 10-6 
QSF star formation efficiency (13.5) 0.07 0.05 — 0.15 
€disk SN feedback disk reheating efficiency ({B.6) 3.5 1-5 
€halo SN feedback halo ejection efficiency (B5) 0.35 0.1 — 0.5 
^fej ejected gas reincorporation efficiency ($3.6) 0.5 0.1 — 1.0 
Tmerger major merger mass ratio threshold (Ba 0.3 0.2 — 0.4 
R instantaneous recycled fraction of SF to the cold disk ($3.9) 0.3 0.2 — 0.4 
Y yield of metals produced per unit SF (BJ 0.03 0.02 — 0.04 


In our model we aim to motivate each aspect of the 
physics of galaxy formation using the best available obser- 
vations and simulations. Our parameters have been chosen 
to reproduce local galaxy properties and are stable in the 
sense that none of our results is critically dependent on any 
single parameter choice; plausible changes in one parame- 
ter or recipe can usually be accommodated through adjust- 
ment of the remaining parameters within their own plausible 
range. The particular model we present is thus not unique. 
Importantly, our model for radio galaxy heating in cooling 
flows, which is the main focus of this paper, is only weakly 
influenced by the remaining galaxy formation and black hole 
growth physics. This is because our radio mode feedback is 
active only in massive objects and at late times and it has 
no effect during the principal growth phase of most galaxies 
and AGN. 

The distribution of galaxy light in our ‘best’ model is 
shown in the bottom panel of Fig. [I] for comparison with 
the mass distribution in the top panel. For both the vol- 
ume is a projected 330 x 280 x 15 !Mpc slice cut from 
the full 0.125 h^? Gpc? simulation box. The plot of surface 
brightness is colour-coded by the luminosity-weighted mean 
B—V colour of the galaxies. On large scales light clearly fol- 
lows mass, but non-trivial biases become evident on smaller 
scales, especially in 'void' regions. The redder colour of 
galaxies in high density regions is also very clear. 


3.2 Gas infall and cooling 


We follow the standard paradigm set out by[White & Fren 
as adapted for implementation on high resolu- 
tion N-body simulations by and 
(2004). This assumes that as each dark mat- 
ter halo collapses, its own ‘fair share’ of cosmic baryons 
collapse with it (but see Section [3.3] below). Thus in our 
model the mass fraction in baryons associated with every 
halo is taken to be fy = 17%, consistent with the first-year 
WMAP result pre nat. Initially these baryons 
are in the form of diffuse gas with primordial composition, 
but later they include gas in several phases as well as stars 
and heavy elements. The fate of the infalling gas depends on 


redshift and on the depth of the halo potential (Silki/1977; 
At late times and in massive systems the gas shocks to the 
virial temperature and is added to a quasi-static hot atmo- 
sphere that extends approximately to the virial radius of the 
dark halo. Gas from the central region of this atmosphere 
may accrete onto a central object through a cooling flow. 
At early times and in lower mass systems the infalling gas 
still shocks to the virial temperature but its post-shock cool- 
ing time is sufficiently short that a quasi-static atmosphere 
cannot form. Rather the shock occurs at much smaller ra- 
dius and the shocked gas cools rapidly and settles onto a 
central object, which we assume to be a cold gas disk. This 
may in turn be subject to gravitational instability, leading 
to episodes of star formation. 

More specifically, the cooling time of a gas is conven- 
tionally taken as the ratio of its specific thermal energy to 
the cooling rate per unit volume, 


3. fim,kT 
teso =en A A . 
! CY p(y MTZ) (3) 


Here pm, is the mean particle mass, k is the Boltzmann 
constant, pg(r) is the hot gas density, and A(T, Z) is the 
cooling function. The latter depends both on the metallicity 
Z and the temperature of the gas. In our models we assume 
the post-shock temperature of the infalling gas to be the 
virial temperature of the halo, T = 35.9 (VA /km s !yK. 
When needed, we assume that the hot gas within a static 
atmosphere has a simple ‘isothermal’ distribution, 
Mhot 


pa(r) = dra , (4) 


where mhot is the total hot gas mass associated with the 
halo and is assumed to extend to its virial radius Ryir. 

'To estimate an instantaneous cooling rate onto the cen- 
tral object of a halo, given its current hot gas content, we 
define the cooling radius, rcoo41, as the radius at which the lo- 
cal cooling time (assuming the structure of Eq.[4) is equal to 
a suitably defined age for the halo. (1991) 
took this age to be the Hubble time ty, while 


(1994 ) used the time scale over which the main pro- 
genitor last doubled its mass. [Somerville & Primacki (1999) 
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systems in the systems that have 


"rapid cooling" formed a static 


infall regime 


hot halo 


a 


Myir / Mo 


Figure 2. The ratio of the cooling radius to virial radius for a random selection of virialised systems at z —0 plotted against their dark 
matter virial mass. Systems identified to be in the ‘rapid cooling’ regime are shown in the left panel, while those that have formed a 
static hot halo are shown on the right (Section[8.2). A sharp transition between the two regimes is seen close to that found by 


(2004), marked by the solid vertical line. 


argued that the time interval since the last major merger 
is more appropriate since such mergers redistribute hot gas 
within the halo. Here we slow Soemeel et all and 
and define the cooling radius as the 
point where the local cooling time is equal to the halo dy- 
namical time, Ryir/Vvir = 0.1 H(z) !. This is an order of 
magnitude smaller than ty and so results in substantially 
smaller cooling radii and cooling rates (typically by a factor 
of 3) than the assumption of (1991). Our 
choice is justified by the tests of (2002) who 
verified explicitly that it results in good object-by-object 
agreement between the amount of gas predicted to con- 
dense in galaxy mass haloes and the amount which actually 
condensed in their high-resolution SPH simulations of the 
formation of a galaxy cluster and its environment. These 
tests assumed primordial abundances in the cooling func- 
tion. When we implement a chemical enrichment model con- 
sistent with the observed element abundances in intracluster 
gas (see Section [3.9] cooling rates in galaxy mass haloes are 
substantially enhanced and we find (as did 
pona that a smaller coefficient than used in the original 
cooling model is required to avoid excessive 
condensation of gas. 

Using the above definition, a cooling rate can now be 
determined through a simple continuity equation, 


Meool = Ar pg (Tcool)TAoolfcool . (5) 


Despite its simplicity, this equation is a good approximation 
to the rate at which gas is deposited at the centre in the 
similarity solution for a cooling flow. 
Putting it all together we take the cooling rate within a 
halo containing a hot gas atmosphere to be 


T'cool Voir 


R2 f 


vir 


(6) 


Mcool = 0.5 mnot 


We assume this equation to be valid when reoo. < Rvir. 
This is the criterion which [White & Freni proposed 
to separate the static hot halo regime from the rapid cooling 
regime (see below). 

In low mass haloes or at high redshifts the formal cool- 


ing radius lies outside the virial radius reoot > Rvir. The 
post-shock gas then cools in less than one sound crossing 
time and cannot maintain the pressure needed to support an 
accretion shock at large radius. The high-resolution spheri- 
cal infall simulations Eea hing & Whied show 
that in this situation the accretion shock moves inwards, 
the post-shock temperature increases and the mass stored in 
the post-shock hot atmosphere decreases, because the post- 
shock gas rapidly cools onto the central object. In effect, all 
infalling material is accreted immediately onto the central 
disk. In this rapid cooling regime we therefore set the cooling 
rate onto the central object to be equal to the rate at which 
new diffuse gas is added to the halo. 


3.2.1 Rapid cooling or cold accretion? 


Although much simplified, the above model was shown by 
ond [Bil et al] 2003) to give reason- 
ably accurate, object-by-object predictions for the cooling 
and accumulation of gas within the galaxies that formed 
in their N-body+SPH simulations. These neglected star- 
formation and feedback effects in order to test the cooling 
model alone. They also assumed primordial abundances in 
the cooling function. In Fig. P] we show the ratio rcoo1/ Ryir 
as a function of virial mass for haloes in the ‘rapid cool- 
ing’ regime (left panel) and in the ‘static halo’ regime (right 
panel) at z — 0 for our ‘best’ galaxy formation model. The 
two regimes are distinguished by the dominant gas com- 
ponent in each halo: when the mass of hot halo gas ex- 
ceeds that of cold disk gas, we say the galaxy has formed 
a static halo, otherwise the system is taken to be in the 
rapid cooling phase. Many haloes classified as ‘rapidly cool- 
ing’ by this criterion have reoo1 < Rvir, which would appar- 
ently indicate a static hot halo. This is misleading, how- 
ever, as systems where cooling is rapid deposit infalling 
gas onto the central galactic disk on a short timescale 
until they have a low-mass residual halo which satisfies 
Tcoo] ^v Ryir. This then persists until the next infall event. 
Averaging over several cycles of this behaviour, one finds 
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that the bulk of the infalling gas cools rapidly. This is why 
we choose to classify systems by their dominant gas com- 
ponent. Note also that a massive hot halo forms imme- 
diately once cooling becomes inefficient, just_as in the 1- 


D infall simulations of |Forcada-Miré & Whitd (1997) and 
Birnboim & Deke . Our classification is thus quite 


robust. 


The transition between the ‘rapid cooling’ and ‘static 
halo’ regimes is remarkably well defined. At z = 0 it oc- 
curs at a halo virial mass of 2-3 x 1014 Mo, and is approx- 
imately independent of redshift out to at least z= 6. This 
is close to the transition mass found for the same cosmol- 
ogy by [Birnboim & Dera] (ond) using spherically symmet- 
ric simulations, and by using fully 3-D 
simulations. This is, in fact, a coincidence since both sets of 
authors assume cooling functions with no metals, whereas 
our ‘best’ model includes heavy elements which substantially 
enhance cooling at the temperatures relevant for the transi- 
tion. Enrichment in the ‘rapid cooling’ regime arises from the 
mixing of infalling primordial gas with supernova-enriched 
gas ejected in a galactic wind (see SectionB.6). If we modify 
our model to assume a zero-metal cooling function, we find 
the transition mass to shift downwards by about a factor of 


2-3, resulting in a lower cooling rate in comparison to the 
work of [Birnboim & Dekel 200d) and (2003). 

'The reason for the more efficient (we would argue overly 
efficient) cooling appears to be different in these two studies. 
The spherical infall simulations of 
showed good agreement with a transition at the 
point predicted using the original cooling radius definition 
of rather than our revised definition 
which was checked explicitly against SPH simulations by 
(2003). Spherical models thus predict more ef- 
ficient cooling than occurs in typical 3-D situations. This ex- 
plains, perhaps, why find a higher 
transition mass than we would predict for zero metallic- 
ity. also showed that the density es- 
timation in the implementation of SPH used elie al 
leads to overcooling of gas in galaxy mass objects as 
compared to their own entropy and energy conserving SPH 
scheme; the effect is large enough to explain shut erc 
find a higher transition mass than we find for their 
assumed cooling function. 


Both and refer to the 
‘rapid cooling’ regime as ‘cold infall'. This is, in fact, a mis- 
nomer. In this mode the accretion shock occurs closer to the 
central object, and so deeper in its potential well than when 
there is a static hot halo. As a result, the pre-shock veloc- 
ity of infalling gas is greater, resulting in a larger post-shock 
temperature. The two modes do not differ greatly in the tem- 
perature to which infalling gas is shocked, but rather in how 
long (compared to the system crossing time) the gas spends 
at the post-shock temperature before its infall energy is lost 
to radiation. Finally, we note that the existence and impor- 
tance of these two modes was the major insight of the origi- 
val work offi (97; [Bine (197A) and ees & Ostiked 
and has been built into all modern theories for galaxy 
formation. A detailed discussion can be found, for example, 


in [White & Frenki (1991). 


3.3 Reionization 


Accretion and cooling in low mass haloes is required to 
be inefficient to explain why dwarf galaxies contain a rela- 
tively small fraction of all condensed baryons 
EA. This inefficiency may in part result from photoionisa- 
tion heating of the intergalactic medium (IGM) which sup- 
presses_the concentration of baryons in shallow potentials 
[Emtathiodliood). More recent models identify the possible 
low-redshift signature of such heating in the faint end of 
the galaxy luminosity function, most notably in the abun- 


dance of the dwarf satellite galaxies in the local group (e.g. 
Tully et_al||2002: 2002). 
(2000) showed that the effect of photoioniza- 


tion heating on the gas content of a halo of mass Myir can 
be modelled by defining a characteristic mass scale, the so 
called filtering mass, Mr, below which the gas fraction fe is 
reduced relative to the universal value: 


jomit 
b 


halo 
Mai) E 
fe (z, Mus) = (n 526 MeMa 


(7) 
The filtering mass is a function of redshift and changes most 
significantly around the epoch of reionization. It was esti- 
mated by Gnedin using high-resolution SLH-P?M simula- 
tions (but see [Hoeft et ants Kravtsov etal] 2004) pro- 
vided an analytic model for these results which distinguishes 
three ‘phases’ in the evolution of the IGM: z > zo, the epoch 
before the first HII regions overlap; zo < z < Zr, the epoch 
when multiple HII regions overlap; z < zr, the epoch when 
the medium is almost fully reionized. They find that choos- 
ing zo=8 and zr =7 provides the best fit to the numerically 
determined filtering mass. We adopt these parameters and 
keep them fixed throughout our paper. This choice results in 
a filtering mass of 4 x 10°Mo at z= zr, and 3 x 10'°Mo b 

the present day. See Appendix B of [Kravtsov et all 


for a full derivation and description of the analytic model. 


3.4 Black hole growth, AGN outflows, and 
cooling suppression 


There is a growing body of evidence that active galactic 
nuclei (AGN) are a critical piece in the galaxy formation 
puzzle. Our principal interest in this paper is their role in 
suppressing cooling flows, thereby modifying the luminosi- 
ties, colours, stellar masses and clustering of the galaxies 
that populate the bright end of the galaxy luminosity func- 
tion. To treat this problem, we first need a physical model 
for the growth of black holes within our galaxies. 


3.4.1 The ‘quasar mode’ 


In our model (which is based closely on that of 
supermassive black holes 
grow during galaxy mergers both by merging with each other 
and by accretion of cold disk gas. For simplicity, black hole 
coalescence is modelled as a direct sum of the progenitor 
masses and thus ignores gravitational wave losses (including 
such losses is found to have little effect on the properties of 
the final galaxy population). We assume that the gas mass 
accreted during a merger is proportional to the total cold 
gas mass present, but with an efficiency which is lower for 
smaller mass systems and for unequal mergers. Specifically, 
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fou Meold (8) 
1 + (280kms~!/Vyir)? ' 


where we have changed the original parameterisation by tak- 
ing 


Amsu,qQ = 


fot = feu (Msat / central) . (9) 


Here fau ~ 0.03 is a constant and is chosen to reproduce the 
observed local MBH —Mbpulge relation (Magorrian et al. 


. In contrast to 
Kauffmann & Haehnelt we allow black hole accre- 


tion during both major and minor mergers although the ef- 
ficiency in the latter is lower because of the mia: /Mcentral 
term. Thus, any merger-induced perturbation to the gas 
disk (which might come from a bar instability or a merger- 
induced starburst — see Section [3.7] can drive gas onto the 
central black hole. In this way, minor merger growth of the 
black hole parallels minor merger growth of the bulge. The 
fractional contribution of minor mergers to both is typically 
quite small, so that accretion driven by major mergers is the 
dominant mode of black hole growth in our model. We refer 
to this as the ‘quasar mode’. [Note that a more schematic 
treatment of black hole growth would suffice for the pur- 
poses of this paper, but in and in 
future work we wish to examine the build-up of the black 
hole population within galaxies in considerably more detail.) 

There is substantial evidence for strong hydrodynamic 


and radiative feedback from optical/UV and X-ray AGN 
(Arav et alJ {2001}; [de Kool et al] {2001}; 
2003). We have not yet explicitly incor- 


porated such feedback in our modelling, and it may well 
turn out to be important (see, for example, the recent 


simulations of |Di Matteo et al/|2005; |Springel et al/|2005al; 
Hopkins et al . We assume ‘quasar mode’ accretion 


to be closely associated with starbursts, so this feedback 
channel may be partially represented in our models by an 
enhanced effective feedback efficiency associated with star 
formation and supernovae (see Section B.6). 


3.4.2 The ‘radio mode’ 


In our model, low energy ‘radio’ activity is the result of 
hot gas accretion onto a central supermassive black hole 
once a static hot halo has formed around the black hole's 
host galaxy. We assume this accretion to be continual and 
quiescent and to be described by a simple phenomenological 
model: 


' MBH hot Wir 

le sen ( sag) (o5 Jemen) NES 
where mpu is the black hole mass, fhot is the fraction of 
the total halo mass in the form of hot gas, Wir cx p is 
the virial velocity of the halo, and kaan is a free parameter 
with units of Mayr ! with which we control the efficiency 
of accretion. We find below that Kaan = 6 x 10795Msyr-! 
accurately reproduces the turnover at the bright end of the 
galaxy luminosity function. Note that fuo: V3,tu is propor- 
tional to the total mass of hot gas, so that our formula is 
simply the product of the hot gas and black hole masses mul- 
tiplied by a constant efficiency and divided by the Hubble 
time ty. In fact, we find fhot to be approximately constant 
for Vir & 150kms~!, so the dependence of mMpu,R on this 
quantity has little effect. The accretion rate given by Eq. [0] 
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Figure 3. The black hole accretion rate density, mpy, as a 
function of redshift for both the ‘quasar’ and the ‘radio’ modes 
discussed in Section [8.4] This figure shows that the growth of 
black holes is dominated by the ‘quasar mode’ at high redshift and 
falls off sharply at z < 2. In contrast, the ‘radio mode’ becomes 
important at low redshifts where it suppresses cooling flows, but 
is not a significant contributor to the overall black hole mass 
budget. 


is typically orders-of-magnitude below the Eddington limit. 
In SectionB]we discuss physical accretion models which may 
reproduce this phenomenology. 

We assume that ‘radio mode’ feedback injects sufficient 
energy into the surrounding medium to reduce or even stop 
the cooling flow described in Section B.2] We take the me- 
chanical heating generated by the black hole accretion of 


Eq. [I0]to be 
Lau — nnnc , (11) 


where 7 = 0.1 is the standard efficiency with which mass 
is assumed to produce energy near the event horizon, and c 
is the speed of light. This injection of energy compensates 
in part for the cooling, giving rise to a modified infall rate 


(Eq. [6) of 


Lan 
n2 


vir 


(12) 


Ri á 
Mecool = Mcool — 


For consistency we never allow m/o) to fall below zero. It is 
worth noting that moo X pu A(Vax)!2 Và. jl (Eq. [6] 
and Theat = 2Lgu/ V4, e mnn fnot Voir (Eq. [I2], so that 


i 1/2 
Mheat MBH ty (13) 
5 1/2 E 
Mcool an A(Voir)4/2 Voir 


(These scalings are exact for an EdS universe; we have omit- 
ted weak coefficient variations in other cosmologies.) Thus 
in our model the effectiveness of radio AGN in suppressing 
cooling flows is greatest at late times and for large values 
of black hole mass. This turns out to be the qualitative be- 
haviour needed for the suppression of cooling flows to repro- 
duce successfully the luminosities, colours and clustering of 
low redshift bright galaxies. 
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Figure 4. The black hole-bulge mass relation for model galaxies 
at the present day. The local observational result in Gl 
is given by the solid line, where the dashed box shows the 
approximate range over which their fit was obtained. 


3.4.9 The growth of supermassive black holes 


Fig. B]shows the evolution of the mean black hole accretion 
rate per unit volume averaged over the entire Millennium 
Simulation. We separate the accretion into the ‘quasar’ and 
‘radio’ modes described above (solid and dashed lines re- 
spectively). Black hole mass growth in our model is dom- 
inated by the merger-driven ‘quasar mode’, which is most 
efficient at redshifts of two to four, dropping by a factor 
of five by redshift zero. This behaviour has similar form to 
but is weaker than the observed evolution with redshift of 


the bright quasar population (e.g. 1990). 
(See also the discussion in 2000). In 
contrast, our ‘radio mode’ is significant only at late times, as 
expected from the scaling discussed above, and for the high 
feedback efficiency assumed in Eq.[IT]it contributes only 5% 
of the final black hole mass density. We will show, however, 
that the outflows generated by this accretion can have a 
major impact on the final galaxy properties. Finally, inte- 
grating the accretion rate density over time gives a present 
day black hole mass density of 3 x 10? M; Mpc ?, consistent 
with recent observational estimates 
[Merlonil[2004). 

The relationship between black hole mass and bulge 
mass is plotted in Fig. A] for the local galaxy population in 
our ‘best’ model. In this figure, the solid line shows the best 
fit to the observations given by [Haring & Rix (2004) for a 
sample of 30 nearby galaxies with well measured bulge and 
black hole masses. Their results only probe masses over the 
range bounded by the dashed lines. Our model galaxies pro- 
duces a good match to these observations with comparable 
scatter in the observed range (see their Fig. 2). 


3.5 Star formation 


We use a simple model for star formation similar to those 
adopted by earlier authors. We assume that all star for- 
mation occurs in cold disk gas, either quiescently or in a 


burst (see Section [3.7]. Based on the observational work of 
(1998), we adopt a threshold surface density for 
the cold gas below which no stars form, but above which gas 
starts to collapse and form stars. According to Keun 
(1996), this critical surface density at a distance R from the 
galaxy centre may be approximated by 


Wir R E —2 
——— || — M 14 
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We convert this critical surface density into a critical mass 
by assuming the cold gas mass to be evenly distributed over 
the disk. The resulting critical cold gas mass is: 


m 9 Wir i Tdisk ) 
men = 38x 10 (ss) (pee) Mo» 09 


where we assume the disk scale length to be rs = (A/V2) Rvir 
(Mo et al.11998), and set the outer disk radius to raisk = 31s 


based on the properties of the Milky Way 
[2000). Here A is the spin parameter of the dark halo in which 
the galaxy resides [secs et adl bong), measured directly 
from the simulation at each time step. When the mass of 


cold gas in a galaxy is greater than this critical value we 
assume the star formation rate to be 


Sais (R) = 120 ( 


Ms = agr (Meola — Merit) / tayn,disk ; (16) 


where the efficiency agr is typically set so that 5 to 15 per- 
cent of the gas is converted into stars in a disk dynamical 
time tayn,disk, which we define to be raisk/Vvir. This star 
formation model produces a global star formation history 
consistent with the observed star formation density of the 
universe out to at least z=2, as shown in Fig. [5] (note that 
this figure also includes star formation through starbursts — 
see Section B7 . 

When implemented in our model, Eq. leads to 
episodic star formation that self-regulates so as to main- 
tain the critical surface density of Eq. [T4] This naturally 
reproduces the observed spiral galaxy gas fractions with- 
out the need for additional parameterisation, as we demon- 
strate in the top panel of Fig. [o] using model Sb/c galaxies 
identified as objects with bulge-to-total luminosity: 1.5 « 
Mi buige — Mrotai < 2.5 (bulge formation is described in 


Section B.7). 


3.6 Supernova feedback 


As star formation proceeds, newly formed massive stars 
rapidly complete their evolution and end their life as super- 
novae. These events inject gas, metals and energy into the 
surrounding medium, reheating cold disk gas and possibly 
ejecting gas even from the surrounding halo. 

The observations of suggest modelling 
the amount of cold gas reheated by supernovae as 


Amreheated = €disk AM ; (17) 


where Am. is the mass of stars formed over some finite time 
interval and eaisk is a parameter which we fix at €disk = 3.5 
based on the observational data. The energy released in this 
interval can be approximated by 
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Figure 5. The star formation rate density of the universe as a 
function of redshift. The symbols show a compilation of observa- 
tional results taken from Fig. 12 PH (20039). 
The solid line shows our ‘best’ model, which predicts that galax- 
ies form much of their mass relatively early. The dashed line (and 


right axis) indicate the increase in stellar mass with redshift. Ap- 
proximately 5096 of all stars form by z —2. 


AEsn = 0.5 chalo Am«Vén , (18) 


where 0.5 Vgy is the mean energy in supernova ejecta per 
unit mass of stars formed, and €naio parametrises the effi- 
ciency with which this energy is able to reheat disk gas. 
Based on a standard initial stellar mass function and stan- 
dard supernova theory we take Vsn = 630 kms^!. In ad- 
dition, for our ‘best’ model we adopt e€naio = 0.35. If the 
reheated gas were added to the hot halo without changing 
its specific energy, its total thermal energy would change by 


AF hot = 0.5 AMreheated Var $ (19) 


Thus the excess energy in the hot halo after reheating is just 
AE excess = AEsn — AEy. When AE excess < 0 the energy 
transferred with the reheated gas is insufficient to eject any 
gas out of the halo and we assume all hot gas remains as- 
sociated with the halo. When excess energy is present, i.e. 
AE excess >0, we assume that some of the hot gas is ejected 
from the halo into an external ‘reservoir’. Specifically, we 
take 


A A Eexcess 
Mejected = ——— — — hot = | €halo 


SN is A * 2 
Ehot ea x) is. ( 0) 


where Enot = 0.5 my V2. is the total thermal energy of the 
hot gas, and we set Amejectea = 0 when this equation gives 
negative values (implying A Fexcess <0 as discussed above). 
This is similar to the traditional semi-analytic feedback 
recipe, Amejectea X Am, / V2, but with a few additions. 
For small Wir the entire hot halo can be ejected and then 
Amejectead must saturate at Amreheatea. Conversely, no hot 
gas can be ejected from the halo for VA. > €halo Ven /€disk; 
ie. for halo circular velocities exceeding about 200kms ! 
for our favoured parameters. 

Ejected gas leaves the galaxy and its current halo in a 
wind or 'super-wind', but it need not be lost permanently. 
As the dark halo grows, some of the surrounding ejecta may 
fall back in and be reincorporated into the cooling cycle. We 
follow Springel et al] 2003) and [De Lucia et all (004) ana 


model this by assuming 
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Figure 6. Selected results for Sb/c galaxies (identified by bulge- 
to-total luminosity, see Section[B.5) for our best model. (top) Gas 
fractions as a function of B-band magnitude. The solid line is a 
representation of the mean behaviour in the (incomplete) sample 
of [Garnett (2003). (middle) The Tully-Fisher relation, where the 
disk circular velocity, Ve, is approximated by Vyir for the dark 
halo. The solid line with surrounding dashed lines represents the 
mean result and scatter found by [Giovanelli et all my d (bot- 
tom) Cold gas metallicity as a function of total stellar mass. The 


solid line represents the result of (2004). 
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Mejected = —Yej Mejected / tayn > (21) 


where Yej controls the amount of reincorporation per dy- 
namical time; typically we take ^e; = 0.3 to 1.0. Such values 
imply that all the ejected gas will return to the hot halo in 
a few halo dynamical times. 

The prescriptions given in this section are simple, as 
well as physically and energetically plausible, but they have 
little detailed justification either from observation or from 
numerical simulation. They allow us to track in a consis- 
tent way the exchange of each halo's baryons between our 
four phases (stars, cold disk gas, hot halo gas, ejecta), but 
should be regarded as a rough attempt to model the complex 
astrophysics of feedback which will surely need significant 
modification as the observational phenomenology of these 
processes is explored in more depth. 

Supernova feedback and star formation act together to 
regulate the stellar growth of a galaxy. This is especially im- 
portant for L « L* galaxies, where feedback can eject most 
of the baryons from the system, reducing the supply of star- 
forming material for time periods much longer than the cool- 
ing/supernova heating cycle. In the middle panel of Fig. [0] 
we plot the Tully-Fisher relation for model Sb/c galaxies 
(see Section [3.5]. The Tully-Fisher relation is strongly in- 
fluenced by the link between star formation and supernova 
heating. The circular velocity of a galactic disk is (to first 
order) proportional to the virial velocity of the host dark 
matter halo and thus to its escape velocity. In our model 
(and most others) this is closely related to the ability of the 
galaxy to blow a wind. The luminosity of a galaxy is deter- 
mined by its ability to turn its associated baryons into stars. 
'The overall efficiency of this process in the face of supernova 
and AGN feedback sets the amplitude of the Tully-Fisher re- 
lation, while the way in which the efficiency varies between 
systems of different circular velocity has a strong influence 
on the slope. 

To predict a Tully-Fisher relation for our model we need 
to assign a maximum rotation velocity to each of our galaxy 
disks. For central galaxies we simply take this velocity to be 
the Vir of the surrounding halo, while for satellite galax- 
ies we take it to be the Vii, of the surrounding halo at the 
last time the galaxy was a central galaxy. This is a crude 
approximation, and for realistic halo structures it is likely 
to be an underestimate both because of the concentration 
of the dark matter distribution and because of the grav- 


itational effects of the baryons (e.g. 1998). Ob- 


taining a good simultaneous fit to observed luminosity func- 
tions and Tully-Fisher relations remains a difficult problem 
within the ACDM paradigm (see, for example 
[2000). Our unrealistic assumption for the disk rotation ve- 
locity actually produces quite a good fit to the observational 
data ofli et all foot, demonstrating that our sim- 
ple star formation and feedback recipes can adequately rep- 
resent the growth of stellar mass across a wide range of 
scales. We find clear deviations from power law behaviour 
for log W < 2.3 (approximately V. £ 100kms~'), where the 
efficiency of removing gas from low mass systems combines 
with our threshold for the onset of star formation to reduce 
the number of stars that can form. The resulting downward 
bend is qualitatively similar to that pointed out in real data 


by |McGaugh et al (2009). These authors show that includ- 


ing the gaseous component to construct a ‘baryonic’ Tully- 


Fisher relation brings the observed points much closer to a 
power-law, and the same is true in the model we present 
here. Limiting star formation in galaxies that inhabit shal- 
low potentials has a strong effect on the faint-end of the 
galaxy luminosity function, as will be seen in Section. [2] 


3.7 Galaxy morphology, merging and starbursts 


In the model we discuss here, the morphology of a galaxy 
is assumed to depend only on its bulge-to-total luminosity 
ratio, which in turn is determined by three distinct physical 
processes: disk growth by accretion, disk buckling to produce 
bulges, and bulge formation through mergers. We treat disk 
instabilities using the simple analytic stability criterion of 


(1998); the stellar disk of a galaxy becomes unsta- 


ble when the following inequality is met, 


Ve 
— l 
(Gmp/rp)!/? 


where we again approximate the rotation velocity of the disk 
Ve by Wir. For each galaxy at each time-step we evaluate the 
left-hand side of Eq. and if it is smaller than unity we 
transfer enough stellar mass from disk to bulge (at fixed rp) 
to restore stability. 

Galaxy mergers shape the evolution of galaxies, affect- 
ing both their morphology and (through induced starbursts) 
their star formation history. Mergers can occur in our model 
between the central galaxy of a dark halo or subhalo and a 
satellite galaxy which has lost its own dark subhalo. Sub- 
structure is followed in the Millennium Run down to a 
20 particle limit, which means that the orbit of a satellite 
galaxy within a larger halo is followed explicitly until its sub- 
halo mass drops below 1.7 x 1019 A^! Mc. After this point, 
the satellite's position and velocity are represented by those 
of the most bound particle of the subhalo at the last time it 
was identified. At the same time, however, we start a merger 


‘clock’ and estimate a merging time for the galaxy using the 

dynamical friction formula of (1987), 
Voir oat 

Gmsat nA ` 


This formula is valid for a satellite of mass Msat orbiting 
in an isothermal potential of circular velocity Vvir at ra- 
dius rasa... We take Msat and Tsat to be the values mea- 
sured for the galaxy at the last time its subhalo could 
be identified. The Coulomb logarithm is approximated by 
In A = In(1 + Mvir/Msat). The satellite is then merged with 
the central galaxy a time tfriction after its own subhalo was 
last identified. If the main halo merges with a larger system 
before this occurs, a new value for tériction is calculated and 
the merger clock is restarted. 

The outcome of the merger will depend on the baryonic 
mass ratio of the two progenitors . When one dominates the 
process, i.e. a small satellite merging with a larger central 
galaxy, the stars of the satellite are added to the bulge of the 
central galaxy and a minor merger starburst (see below) will 
result. The cold gas of the satellite is added to the disk of 
the central galaxy along with any stars that formed during 
the burst. Such an event is called a minor merger. 

If, on the other hand, the masses of the progenitors are 
comparable a major merger will result. Under these circum- 
stances the starburst is more significant, with the merger 


(22) 


friction = 1.17 (23) 
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Figure 7. The mean condensation rate, (fhcool) as a function of halo virial velocity Vyir at redshifts of 6, 3, 1, and 0. Solid and dashed 
lines in each panel represent the condensation rate with and without ‘radio mode’ feedback respectively, while the vertical dotted lines 
show the transition between the rapid cooling and static hot halo regimes, as discussed in Section [3.2] This figure demonstrates that 
cooling flow suppression is most efficient in our model for haloes with Vj, > 150 kms-! and at z <3. 


destroying the disks of both galaxies to form a spheroid in 
which all stars are placed. The dividing line between a ma- 
jor and minor merger is given by the parameter Tmerger: 
when the mass ratio of the merging progenitors is larger 
than Tmerger à major merger results, otherwise the event is 
a minor merger. Following Sorinel ot all (2001al) we choose 
Tmerger = 0.3 and keep this fixed throughout. 

Our starburst implementation is based on the ‘colli- 
sional starburst’ model of [Somerville et al] (2001). In this 


model, a fraction epurst of the combined cold gas from the 
two galaxies is turned into stars as a result of the merger: 


\ Burst ] (24) 


Cburst — Bourst (Msat /Mcentral 


where the two parameters are taken as Qburst = 0.7 and 
Bourst = 0.56. This model provides a good fit to the numeri- 


cal results of [Cox et al] (2004!) and also|Mihos & Hernquist 
(1994 1996) for merger mass ratios ranging from 1:10 to 1:1. 


3.8 Spectroscopic evolution and dust 


The photometric properties of our galaxies are cal- 
culated using stellar population synthesis models from 
(2003). Our implementation is fully de- 
scribed in and we refer the reader 
there (and to references therein) for further details. 

To include the effects of dust when calculating galaxy 
luminosities we apply the simple ‘plane-parallel slab’ model 
of (1999). 'This model is clearly oversim- 
plified, but it permits us to make a reasonable first-order 
correction for dust extinction in actively star-forming galax- 
ies. For the details of this model we refer the reader to 


Kauffmann et al| (1999) and to references therein. 


3.9 Metal enrichment 


Our treatment of metal enrichment is essentially identical 
to that described in|De Lucia et al| (2004). In this model a 


yield Y of heavy elements is returned for each solar mass of 
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Figure 8. Galaxy luminosity functions in the K (left) and bj (right) photometric bands, plotted with and without ‘radio mode’ feedback 
(solid and long dashed lines respectively — see Section [3.2]. Symbols indicate observational results as listed in each panel. As can be seen, 
the inclusion of AGN heating produces a good fit to the data in both colours. Without this heating source our model overpredicts the 
luminosities of massive galaxies by about two magnitudes and fails to reproduce the sharp bright end cut-offs in the observed luminosity 


functions. 


stars formed. These metals are produced primarily in the su- 
pernovae which terminate the evolution of short-lived, mas- 
sive stars. In our model we deposit them directly into the 
cold gas in the disk of the galaxy. (An alternative would 
clearly be to add some fraction of the metals directly to 
the hot halo. Limited experiments suggest that this makes 
little difference to our main results.) We also assume that 
a fraction R of the mass of newly formed stars is recycled 
immediately into the cold gas in the disk, the so called *in- 
stantaneous recycling approximation’ (see [Cole et al]2000). 
For full details on metal enrichment and exchange processes 
in our model see[De Lucia et al] (2004). In the bottom panel 
of Fig. [6] we show the metallicity of cold disk gas for model 
Sb/c galaxies (selected, as before, by bulge-to-total luminos- 
ity, as described in Section[B.5] as a function of total stellar 
mass. For comparison, we show the result of [Tremonti et al] 
for mean HII region abundances in SDSS galaxies. 


4 RESULTS 


In this section we examine how the suppression of cooling 
flows in massive systems affects galaxy properties. As we will 
show, the effects are only important for high mass galaxies. 
Throughout our analysis we use the galaxy formation model 
outlined in the previous sections with the parameter choices 
of Table[I]except where explicitly noted. 


4.1 The suppression of cooling flows 


We begin with Fig. [Z] which shows how our ‘radio mode’ 
heating model modifies gas condensation. We compare mean 
condensation rates with and without the central AGN heat- 
ing source as a function of halo virial velocity (solid and 
dashed lines respectively). Recall that virial velocity pro- 
vides a measure of the equilibrium temperature of the sys- 
tem through Tyir x V2,, as indicated by the scale on the top 
axis. The four panels show the behaviour at four redshifts 
between six and the present day. The vertical dotted line in 
each panel marks haloes for which reoo1 = Rwir, the transi- 
tion between systems that form static hot haloes and those 
where infalling gas cools rapidly onto the central galaxy disk 
(see section B.2]and Fig. 2). This transition moves to haloes 
of lower temperature with time, suggesting a ‘down-sizing’ of 
the characteristic mass of actively star-forming galaxies. At 
lower Vyir gas continues to cool rapidly, while at higher Vii. 
new fuel for star formation must come from cooling flows 
which are affected by ‘radio mode’ heating. 


The effect of ‘radio mode’ feedback is clearly substan- 
tial. Suppression of condensation becomes increasingly effec- 
tive with increasing virial temperature and decreasing red- 
shift. The effects are large for haloes with Vi, Z 150kms ! 
(Tir & 109K) at z $3. Condensation stops almost com- 
pletely between z — 1 and the present in haloes with 
Voir  300kms^! (Tyir > 3 x 10°K). Such systems corre- 
spond to the haloes of groups and clusters which are typ- 
ically observed to host massive elliptical or cD galaxies at 
their centres. Our scheme thus produces results which are 
qualitatively similar to the ad hoc suppression of cooling 
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Figure 9. The B—V colours of model galaxies plotted as a func- 
tion of their stellar mass with (top) and without (bottom) ‘radio 
mode’ feedback (see Section [3.2). A clear bimodality in colour is 
seen in both panels, but without a heating source the most mas- 
sive galaxies are blue rather than red. Only when heating is in- 
cluded are massive galaxies as red as observed. Triangles (red) and 
circles (blue) correspond to early and late morphological types re- 
spectively, as determined by bulge-to-total luminosity ratio (see 
Section 4,2). The thick dashed lines mark the resolution limit to 
which morphology can be reliably determined in the Millennium 
Run. 


flows assumed in previous models of galaxy formation. For 
example, switched off gas condensa- 
tion in all haloes with Vir > 350kms !, whilelEintton all 
stopped condensation when the bulge mass exceeded 
a critical threshold. 


4.2 Galaxy properties with and without AGN 
heating 


The suppression of cooling flows in our model has a dramatic 
effect on the bright end of the galaxy luminosity function. 
In Fig. [8] we present K- and bj-band luminosity functions 
(left and right panels respectively) with and without ‘ra- 
dio mode’ feedback (solid and dashed lines respectively). 
The luminosities of bright galaxies are reduced by up to 
two magnitudes when the feedback is switched on, and this 
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Figure 10. Mean stellar ages of galaxies as a function of stellar 
mass for models with and without 'radio mode' feedback (solid 
and dashed lines respectively). Error bars show the rms scatter 
around the mean for each mass bin. The suppression of cooling 
flows raises the mean age of high-mass galaxies to large values, 
corresponding to high formation redshifts. 


induces a relatively sharp break in the luminosity function 
which matches the observations well. We demonstrate this 
by overplotting K-band data from and 
2003) in the left panel, and b;-band data from 
(2002) in the right panel. In both band-passes 
the model is quite close to the data over the full observed 
range. We comment on some of the remaining discrepancies 
below. 

Our feedback model also has a significant effect on 
bright galaxy colours, as we show in Fig.[9] Here we plot the 
B—V colour distribution as a function of stellar mass, with 
and without the central heating source (top and bottom pan- 
els respectively). In both panels we have colour-coded the 
galaxy population by morphology as estimated from bulge- 
to-total luminosity ratio (split at L»uge/Ltotai = 0.4). Our 
morphological resolution limit is marked by the dashed line 
at a stellar mass of ~ 4 x 10? M; this corresponds approx- 
imately to a halo of 100 particles in the Millennium Run. 
Recall that a galaxy's morphology depends both on its past 
merging history and on the stability of its stellar disk in our 
model. Both mergers and disk instabilities contribute stars 
to the spheroid, as described in Section[B.7] The build-up of 
haloes containing fewer than 100 particles is not followed in 
enough detail to model these processes robustly. 

A number of important features can be seen in Fig. J] 
Of note is the bimodal distribution in galaxy colours, with 
a well-defined red sequence of appropriate slope separated 
cleanly from a broader ‘blue cloud’. It is significant that 
the red sequence is composed predominantly of early-type 
galaxies, while the blue cloud is comprised mostly of disk- 
dominated systems. This aspect of our model suggests that 
that the physical processes which determine morphology 
(i.e. merging, disk instability) are closely related to those 
which control star formation history (i.e. gas supply) and 
thus determine galaxy colour. The red and blue sequences 
both display a strong metallicity gradient from low to high 
mass (c.f. Fig. (6), and it is this which induces a ‘slope’ in 
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Figure 11. The b;j-band galaxy luminosity function split by colour at B— V = 0.8 (Fig.[9) into blue (left panel) and red (right panel) sub- 
populations (solid lines). The dotted lines in each panel repeat the opposite colour luminosity function for reference. Symbols indicate the 
observational results of[Madgwick et al] for early and late-type 2dFGRS galaxies, split according to spectral type. Although our 
model split by colour captures the broad behaviour of the observed type-dependent luminosity functions, there are important differences 


which we discuss in Section [1.2] 


the colour-magnitude relations which agrees well with ob- 
servation (e Bad a aloo. 

By comparing the upper and lower panels in Fig. B] we 
can see how ‘radio mode’ feedback modifies the luminosities, 
colours and morphologies of high mass galaxies. Not surpris- 
ingly, the brightest and most massive galaxies are also the 
reddest and are ellipticals when cooling flows are suppressed, 
whereas they are brighter, more massive, much bluer and 
typically have disks if cooling flows continue to supply new 
material for star formation. AGN heating cuts off the gas 
supply to the disk from the surrounding hot halo, truncating 
star formation and allowing the existing stellar population 
to redden. However, these massive red galaxies do continue 
to grow through merging. This mechanism allows the dom- 
inant cluster galaxies to gain a factor 2 or 3 in mass with- 
out significant star formation, in apparent agreement with 
observation [sracon-Salamanca ef afin). This late-stage 
(i.e. z S 1) hierarchical growth moves objects to higher mass 
without changing their colours. 


It is also interesting to examine the effect of AGN heat- 
ing on the stellar ages of galaxies. In Fig. [0] the solid and 
dashed lines show mean stellar age as a function of stel- 
lar mass for models with and without ‘radio mode’ feed- 
back, while error bars indicate the rms scatter around the 
mean. Substantial differences are seen for galaxies with 
Mstetlar © 10!! Mo: the mean age of the most massive galax- 
ies approaching 12 Gyr when cooling flows are suppressed 
but remaining around 8 Gyr when feedback is switched off. 
Such young ages are clearly inconsistent with the old stel- 


lar populations observed in the majority of massive cluster 
ellipticals. 


The colour bimodality in Fig. Plis so pronounced that it 
is natural to divide our model galaxies into red and blue pop- 
ulations and to study their properties separately. We do this 
by splitting at B—V = 0.8, an arbitrary but natural choice. 
Fig. [shows separate b; luminosity functions for the result- 
ing populations. For comparison we overplot observational 
results from for 2dFGRS galaxies 
split by spectral type. Their luminosity functions are es- 
sentially identical to those of (2003), who split 
the 2dFGRS catalogue by b;-rr colour. It thus can serve to 
indicate the observational expectations for populations of 
different colour. The broad behaviour of the red and blue 
populations is similar in the model and in the 2dFGRS. The 
faint-end of the luminosity function is dominated by late- 
types, whereas the bright end has an excess of early-types. 
The two populations have equal abundance about half to 


one magnitude brighter than Mj, (Norberg et a1/|2002). 


Fig. [I1]also shows some substantial differences between 
model and observations. The red and blue populations differ 
more in the real data than they do in the model. There is a 
tail of very bright blue galaxies in the model, which turn out 
to be objects undergoing strong, merger-induced starbursts. 
'These correspond in abundance, star formation rate and evo- 
lutionary state to the observed population of Ultraluminous 
Infrared Galaxies (ULIRG's) with the important difference 
that almost all the luminosity from young stars in the real 
systems is absorbed by dust and re-emitted in the mid- to 


far-infrared (Sanders & Mirabell1998). Clearly we need bet- 
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ter dust modelling than our simple ‘slab’ model (SectionB-3) 
in order to reproduce the properties of such systems ade- 
quately. If we suppress starbursts in bright galaxy mergers 
we find that the blue tail disappears and the observed be- 
haviour is recovered. A second and substantial discrepancy 
is the apparent overproduction of faint red galaxies in our 
model, as compared to the 2dF measurements (however see 
[Popesso ot all 2005 [Gonzalez ot all[2005). Further work is 
clearly needed to understand the extent and significance of 
this difference. 


5 PHYSICAL MODELS OF AGN FEEDBACK 


Our phenomenological model for ‘radio mode’ feedback (Sec- 
tion[.4.2) is not grounded in any specific model for hot gas 
accretion onto a black hole or for the subsequent generation 
and thermalization of energy through radio outflows. Rather 
it is based on the observed properties of cooling flows and 
their central radio sources, and on the need for a source of 
feedback which can suppress gas condensation onto massive 
galaxies without requiring the formation of new stars. We 
have so far focused on the effects of such feedback without 
discussing how it might be realised. In this section we present 
two physical models which suggest how accretion onto the 
central black hole may lead to activity in a way which could 
justify the parameter scalings we have adopted. 


5.1 Cold cloud accretion 


A simple picture for cooling flow evolution, based on the sim- 
ilarity solution of [Bertachinged for an unperturbed 
halo in isolation, can be summarised as follows. Cooling 
flows develop in any halo where the cooling time of the bulk 
of the hot gas is longer than the age of the system so that a 
static hot halo can form. Such haloes usually have a strong 
central concentration and we approximate their structure by 
a singular isothermal sphere. The inner regions then have a 
local cooling time shorter than the age of the system, and 
the gas they contain radiates its binding energy and flows 
inwards. The flow region is bounded by the cooling radius 
Tcoo! Where the local cooling time is equal to the age of the 
system (see section B.2). This radius increases with time as 
1/2. As Bertschinger showed, the temperature of the gas 
increases by about 2096 as it starts to flow inwards, and its 
density profile flattens to pg c pm. Initially, the flow is 
subsonic and each gas shell sinks stably and isothermally 
in approximate hydrostatic equilibrium. As it sinks, how- 
ever, its inward velocity accelerates because its cooling time 
shrinks more rapidly than the sound travel time across it, 
and at the sonic radius, rsonic, the two become equal. At 
this point the shell goes into free fall, its temperature de- 
creases rapidly and it may fragment as a result of thermal 
instability (Cowie ei aL[1980 [Nulscul 98d [Balbus & Sokel 
[1989). 'The dominant component of the infalling gas is then 
in the form of cold clouds and is no longer self-coupled by 
hydrodynamic forces. Different clouds pursue independent 
orbits, some with pericentres perhaps orders of magnitude 
smaller than rsonic. If these lie within the zone of influence of 
the black hole, reu = Gmpu/ Và. we assume that some of 
the cold gas becomes available for fuelling the radio source; 
otherwise we assume it to be added to the cold gas disk. 


The parameter scalings implied by this picture can be 
estimated as follows. The sound travel time across a shell 
at the cooling radius is shorter than the cooling time by a 
factor ~ Tcool/Rvir. At smaller radii the ratio of cooling time 
to sound travel time decreases as r!/? so that l'sonic/ cool ^v 
(Toot / Ri)? implying Trsonic ~ roo) / Rein If we adopt rpu > 
107? rsonic as the condition for effective fuelling of the radio 
source, we obtain 


MBH > 1074 Mvyir (rco01/Rvir)? (25) 


as the corresponding minimum black hole mass for frag- 
mented clouds to be captured. Under such conditions, only 
a small fraction (~ 0.01%) of the cooling flow mass need 
be accreted to halt the flow. The ratio in parentheses on 
the right-hand side of this equation scales approximately 
as Tcool/Rvir (X (Mnot/Mvir) Pta P VI, so the minimum 
black hole mass scales approximately as (mno: / Myir)?/ as d 
and is almost independent of Vvir. In our model, the growth 
of black holes through mergers and ‘quasar mode’ accretion 
produces a population where mass increases with time and 
with host halo mass. As a result, effective fuelling takes place 
primarily in the more massive haloes and at late times for 
this ‘cold cloud’ prescription. 

To test this particular model we switch off our stan- 
dard phenomenological treatment of ‘radio mode’ feedback 
(section 8.4.2), assuming instead that feedback occurs only 
when Eq. Blis satisfied and that in this case it is sufficient 
to prevent further condensation of gas from the cooling flow. 
All other elements of our galaxy and black hole formation 
model are unchanged. The resulting cooling flow suppression 
is similar to that seen in Fig. [7] and all results presented 
in Section B] and [4] are recovered. An illustration of this is 
given by Fig. [2] where we compare the K-band luminos- 
ity function from this particular model (the dashed line) to 
the observational data (c.f. also Fig. [8). The model works so 
well, of course, because the numerical coefficient in Eq. [25]is 
uncertain and we have taken advantage of this to choose a 
value which puts the break in the luminosity function at the 
observed position. This adjustment plays the role of the effi- 
ciency parameter «Kacan in our standard analysis (see Eq.[I0].. 


5.2 Bondi-Hoyle accretion 


Our second physical model differs from the first in assuming 
that accretion is not from the dominant, cold cloud compo- 
nent which forms within the sonic radius, but rather from a 
subdominant hot component which fills the space between 
these clouds. The clouds themselves are assumed to be lost 
to the star-forming disk. The density profile of the residual 
hot component was estimated by ess Fabi] 
from the condition that the cooling time of each radial shell 
should remain equal to the sound travel time across it as it 
flows inwards. This requires the density of the hot compo- 
nent to vary as 1/r within rsonic and thermal instabilities 
must continually convert material into condensed clouds in 
order to maintain this structure as the hot gas flows in. 
'The rate at which hot gas is accreted onto the black 
hole can then be estimated from the Bondi-Hoyle formula 


(Bonai[1952. [Edgar 2004): 


2 
ThiBondi = 2.51 G? DRE. i (26) 


Cs 
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Figure 12. The observed K-band galaxy luminosity function is 
compared with the results from models using the two physical pre- 
scriptions for ‘radio mode’ accretion discussed in Section [B] the 
Bondi-Hoyle accretion model (solid line) and the cold cloud ac- 
cretion model (dashed line). Symbols indicate observational data 


from|Cole et al} and (2003). Both models can 


produce a luminosity function which matches observation well. 


Here po is the (assumed uniform) density of hot gas around 
the black hole, and in all that follows we approximate the 
sound speed, cs, by the virial velocity of the halo, Vir. Of 
course, the density distribution of gas surrounding the black 
hole is not uniform so the question immediately arises as 
to what density we should choose. We follow a suggestion 
of E. Churazov and use the value predicted by the 'maxi- 
mum cooling flow’ model of (2000) at the 
Bondi radius r Bondi = 2G-Mpu/c2 = 2rpu, the conventional 
boundary of the sphere of influence of the black hole. We 
therefore equate the sound travel time across a shell at this 
radius to the local cooling time there: 


2rpondi , 4Gmen _ 3 fiomypkT (27) 


Cs Va. 7 2 Pg (rBenai) A(T, Z) f 
Solving for the density gives 
3UMp kT Va. 
8G A MBH i 
Combining Eq. 28] with [26] provides us with the desired es- 
timate for the hot gas accretion rate onto the black hole: 


po = Pg(TBonai) = (28) 


kT 
Bondi £2 GUMp—— MBH . (29) 


A 


Notice that this rate depends only on the black hole mass 
and on the virial temperature of the halo. It is independent 
both of time and of mnot/Mvir, the hot gas fraction of the 
halo. It is valid as long as rBonai < rsonic, which is always 
the case in our models. 

To investigate the effects of this model we replace the 
phenomenological ‘radio mode’ accretion rate of Eq. [TO] with 
that given by Eq.[29] Since the latter has no adjustable effi- 
ciency, we use the energy generation parameter 7 of Eq. [IT] 
to control the effectiveness of cooling flow suppression. (This 
was not necessary before since 7 always appeared in the 


product 7 &AGN, where Kaan is the efficiency parameter of 
Eq.[f0]) With this change of cooling flow accretion model 
and taking 7 = 0.03, we are able to recover the results 
of Sections B] and B] without changing any other aspects 
of our galaxy and black hole formation model. The final 
galaxy population is, in fact, almost identical to that pre- 
sented in previous sections. This is not surprising, perhaps, 
since Eq. [29]has very similar scaling properties to Eq. [10] In 
Fig.[I2]we illustrate the success of the model by overplotting 
its prediction for the K-band luminosity function (the solid 
line) on the observational data and on the prediction of the 
cold cloud accretion model of the last subsection. The two 
models agree very closely both with each other and with our 
standard phenomenological model (see Fig. Bp. 


6 CONCLUSIONS 


AGN feedback is an important but relatively little explored 
element in the co-evolution of galaxies and the supermassive 
black holes at their centres. In this paper we set up machin- 
ery to study this co-evolution in unprecedented detail using 
the very large Millennium Run, a high-resolution simula- 
tion of the growth of structure in a representative region of 
the concordance ACDM cosmology. Most of our modelling 
follows earlier work, but in an important extension we in- 
troduce a ‘radio’ feedback mode, based on simple physical 
models and on the observed phenomenology of radio sources 
in cooling flows. This mode suppresses gas condensation at 
the centres of massive haloes without requiring the forma- 
tion of new stars. Our modelling produces large catalogues 
of galaxies and supermassive black holes which can be used 
to address a very wide range of issues concerning the evo- 
lution and clustering of galaxies and AGN. Some clustering 
results were already presented in{Springel et al] (2005b). In 
the present paper, however, we limit ourselves to present- 
ing the model in some detail and to investigating the quite 
dramatic effects which ‘radio mode’ feedback can have on 
the luminosities and colours of massive galaxies. Our main 
results can be summarised as follows: 


(i) We study the amount of gas supplied to galaxies in each 
of the two gas infall modes discussed by 
(1991): the ‘static halo’ mode where postshock cooling is 
slow and a quasistatic hot atmosphere forms behind the 
accretion shock; and the ‘rapid cooling’ mode where the 
accretion shock is radiative and no such atmosphere is 
present. We distinguish these two modes using the crite- 
rion of as modified bv emnes et all 
2001a) and tested explicitly using SPH simulations by 
(2002). Our results show a sharp transition 
between the two regimes at a halo mass of 2-3 x 10!! Mọ. 
'This division depends on the chemical enrichment prescrip- 
tion adopted and moves from higher to lower Vvir with time 
(corresponding to approximately constant Myir), suggesting 
a ‘down-sizing’ of star formation activity as the bulk of the 
gas accreted by the haloes of larger systems is no longer 
available in the interstellar medium of the central galaxy. 


(ii) We have built a detailed model for cooling, star for- 
mation, supernova feedback, galaxy mergers and metal en- 
richment based on the earlier models of 
(1999), [Springel et al] and (2004). 
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Applied to the Millennium Run this model reproduces many 
of the observed properties of the local galaxy population: 
the Tully-Fisher, cold gas fraction/stellar mass and cold gas 
metallicity/stellar mass relations for Sb/c spirals (Fig. [6]; 
the field galaxy luminosity functions (Fig. [S] & [TT]; the 
colour-magnitude distribution of galaxies (Fig. 9); and the 
increase in mean stellar age with galaxy mass (Fig. [I0]. In 
addition the model produces a global star formation his- 
tory in reasonable agreement with observation (Fig. 5). We 
also show in |Springel et all that the z = 0 cluster- 
ing properties of this population are in good agreement with 
observations. 


(iii) Our black hole implementation extends the previous 
work of Kaun © Heche by assuming three 
modes of AGN growth: merger-driven accretion of cold disk 
gas in a ‘quasar mode’, merging between black holes, and 
‘radio mode’ accretion when a massive black hole finds itself 
at the centre of a static hot gas halo. The ‘quasar mode’ is 
the dominant source for new black hole mass and is most 
active between redshifts of four and two. The ‘radio mode’ 
grows in overall importance until z=0 and is responsible for 
the feedback which shuts off the gas supply in cooling flows. 
This model reproduces the black hole mass/bulge mass re- 
lation observed in local galaxies (Fig. Æ). The global history 
of accretion in the ‘quasar mode’ is qualitatively consistent 
with the evolution inferred from the optical AGN population 


(Fig. Bp. 


(iv) Although the overall accretion rate is low, ‘radio mode’ 
outflows can efficiently suppress condensation in massive 
systems (Fig.[7). As noted by many authors who have stud- 
ied the problem in more detail than we do, this provides an 
energetically feasible solution to the long-standing cooling 
flow ‘problem’. Our analysis shows that the resulting sup- 
pression of gas condensation and star formation can produce 
luminosity functions with very similar bright end cut-offs to 
those observed (Fig. [S], as well as colour-magnitude dis- 
tributions in which the most massive galaxies are red, old 
and elliptical, rather than blue, young and disk-dominated 


(Figs PB] and [10]. 


(v) The B—V colour distribution of galaxies is bimodal at 
all galaxy masses. Galaxies with early-type bulge-to-disk ra- 
tios are confined to the red sequence, as are the most mas- 
sive galaxies, and the most massive galaxies are almost all 
bulge-dominated, as observed in the real universe (Fig. J). 
'This bimodality provides a natural division of model galax- 
ies into red and blue subpopulations. The colour-dependent 
luminosity functions are qualitatively similar to those found 
for early and late-type galaxies in the 2dFGRS (Fig. [TT], 
although there are significant discrepancies. After exhaust- 
ing their cold gas, massive central galaxies grow on the 
red sequence through ‘burstless’ merging, gaining a factor 
of two or three in mass without significant star formation 
(Aragon-Salamanca et al] [1998). Such hierarchical growth 
does not change a galaxy's colour significantly, moving it 
brightward almost parallel to the colour-magnitude relation. 


(vi) We present two physical models for black hole accretion 
from cooling flow atmospheres. We suppose that this accre- 
tion is responsible for powering the radio outflows seen at 


the centre of almost all real cooling flows. The models differ 
in their assumptions about how gas accretes from the inner 
regions of the cooling flow, where it is thermally unstable 
and dynamically collapsing. One assumes accretion of cold 
gas clouds if these come within the sphere of influence of 
the black hole, while the other assumes Bondi-like accretion 
from the residual diffuse hot gas component. Each of the two 
models can produce z=0 galaxy populations similar both to 
that of our simple phenomenological model for ‘radio mode’ 
feedback and to the observed population (see Fig.[12). Our 
main results are thus not sensitive to the details of the as- 
sumed accretion models. 


The presence of heating from a central AGN has long 
been suspected as the explanation for the apparent lack of 
gas condensation in cluster cooling flows. We have shown 
that including a simple treatment of this process in galaxy 
formation models not only ‘solves’ the cooling flow problem, 
but also dramatically affects the properties of massive galax- 
ies, inducing a cut-off similar to that observed at the bright 
end of the galaxy luminosity function, and bringing colours, 
morphologies and stellar ages into much better agreement 
with observation than is the case for models without such 
feedback. We will extend the work presented here in a com- 
panion paper, where we investigate the growth of supermas- 
sive black holes and the related AGN activity as a function 
of host galaxy properties out to high redshift. The cata- 
logues of galaxies and supermasive black holes produced by 
our modelling machinery are also being used for a very wide 
range of projects related to understanding formation, evo- 
lution and clustering processes, as well as for interpreting 
observational samples. 
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